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^ : ABSTRACT 
O 

We model the thermal effect of young stars on their surrounding environment 
Q I in order to understand clustered star formation. We take radiative heating of 

dust, dust-gas coUisional heating, cosmic-ray heating, and molecular cooling into 
CN I account. Using Dusty, a spherical continuum radiative transfer code, we model 

the dust temperature distribution around young stellar objects with various lu- 
minosities and surrounding gas and dust density distributions. We have created 
a grid of dust temperature models, based on our modeling with Dusty, which we 
can use to calculate the dust temperature in a field of stars with various param- 
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^ I eters. We then determine the gas temperature assuming energy balance. Our 

models can be used to make large-scale simulations of clustered star formation 
more realistic. 
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1. Introduction 



^ ■ Most of the stars in our galaxy form in groups or clusters ( Lada fc Ladall2003 ). There- 



c3 ■ fore, in order to understand the star formation history, the shape of the mass function, and 
the formation of massive (M > 5 Mq) stars in our galaxy, the star formation process must 
be studied in its most common environment - a cluster. As stars form from their initial 
reservoir of gas and dust, they interact with their environment and heat the surrounding 
material, thus affecting future star formation. One of the first effects a protostar has on its 
environment is radiative heating from the accretion luminosity and, subsequently, nuclear 
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fusion. The radiation efficiently heats the dust, which in turn heats the gas through colh- 
sions. Young stars also affect their environment via strong winds and ionization, but this 
only occurs when they are very massive and have evolved past the very early stages of star 
formation. We assume that the massive stars in our sample are very young and are accreting 
at very high rates (M > 1 x lO'^Mfn/yr ). This high accre tion rate allows the infalling mass 
to absorb all of the stellar UV photons (IChurchwelll |2002| ) . 



Many groups use large scale computer simulations to model clustered star formation. 
This is a complicated process requir i ng m any assumptions in order to make the problem 
tractable. iKlessen. Burkert. &: Bate! (119981) and Marte l, Eva. ns, & Shapiro (2006) assume 



that the gas is isothermal. iBate. Bonnell. fc Bromnj (120031 ) go beyond this assumption 



by using a barotropic equation of state. However, until recently, no one has included 
the effect of radiatively heating t he dust and gas by the stars formed in the simulation. 
Krumholz. Klein, fc McKed (120071 ) have included an approximate radiative transfer method, 
which works well in optically thick regions. Their method assumes that the gas temperature 
is equal to the dust temperature throughout their simulation. This approximation is only 
valid at high densities when the dust and gas are coUisionally coupled. Here we develop a 
method that explores the effect of radiative heating and the dust and gas energetics for a 
range of optical depths and densities. 

In our method we include various heating and cooling processes to calculate the dust 
and gas temperature. Stars can heat dust grains more effectively than the gas because dust 
grains have broad-band absorption properties. Although we will not be explicitly modeling 
the motion or energy density of dust grains, we assume the dust and gas are well-mixed and 
the dust grains t ransfer energy to ga s particles through collisions using the energy transfer 
rate discussed in lYoung et al.l (120041 ). The gas is heated by collisions with hot dust grains 
and cosmic rays. It can cool through CO and other molecular line emission. 

In this paper, we calculate the dust and gas temperature in a field of stars. The dust 
temperature around a single source is calculated using a look-up table which we develop 
here. With this look-up table and an approximation to the flux-temperature conversion, we 
calculate the dust temperature in the field. Our look-up table is needed since the calculation 
of a single dust temperature distribution can take longer than a minute on current desktops 
and would take a substantial fraction of a large-scale simulation's computations. Therefore 
we outline our method here which can be used to decrease the time spent on the calculation 
of the dust temperature in future studies of clustered star formation. 

With the calculated value of the dust temperature, we derive the gas temperature field 
for a distribution of stellar sources, as in a young stellar cluster. The effect that protostars 
have on heating their environment using a hydrodynamic and gravity simulation will be 
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addressed in a future paper. In this paper, we first discuss the calculation of the dust 
temperature for single and multiple sources (§2]), then we describe our gas temperature 
calculation (§3]), and finally, we show some dust and gas temperature distributions when 
multiple sources are present (§!]). 



We consider two methods of calculating the dust temperature when there are multiple 
heating sources. The first approach assumes radiative equilibrium after summing up the 
fiux of multiple sources heating a dust grain. This approach makes simplifying assumptions 
about the dust absorption and emission properties. The second approach uses the one- 
dimensional spherical radiative transfer code Dusty (Nenkova, Ivezic, & Elitzur 2000), to add 
up the energy density contributed by multiple sources; the dust temperature distribution is 
calculated separately for each source, then the temperatures are converted to energy densities 
(assuming radiative equilibrium), which are summed together and then converted back to 
a temperature. We first discuss the analytic approach ( §2.ip . then the numerical approach 
to calculating the dust temperature ( §2.21) . and finally we compare and analyze the two 
approaches ( §2.3p . 



where R^i is the radius of star i, S^i is the fiux density at the stellar surface of star i (which 
we assume is a blackbody), Qa is the dust grain's absorption efficiency as a function of 
frequency, vra^ is the projected surface area of the grain exposed to the star's light, and the 
separation between the star i and the dust grain is 

where r^,j is the location of star i and r is the position of the dust grain. Substituting in the 
Planck function for star i at temperature T^i and assuming 



2. Dust Temperature Calculation 



2.1. Analytic Dust Temperature Calculation 



The rate of energy absorbed by a dust grain in a field of stars is 
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equation ([T]) becomes 
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where 



/4+^„ = r(4 + /3,)C(4 + /3a) (4) 
and the functions T{x) and are defined as the gamma and Riemann zeta function. 

The same assumptions can be made about the emission of the grain except that the 
grain emits in all directions (i.e. no? becomes Ana"^), the grain is emitting instead of the star 
(i.e. T^i becomes T^), and the dust grain's emission efficiency is 
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These assumptions give 
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Assuming radiative equilibrium and removing the dependence on stellar radius with L^, 
ATiRlaT^, Trf at position r is given by 
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If we assume that grains absorb like graybodies (/?„ = 0) where stars emit most of their 
luminosity (UV), then Qa{i^) = Qa{UV). The value of Qe is given at 125/im with respect to 
Qa- Using this approximation, we obtain 
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Notice that the dependence on stellar temperature disappears in this case. 

The value of (5a(f^V")/(5e(125/im), where the UV wavelength range is 0.15/im — 0.30/im, 
can be calculated from observations or from a dust grain model. We compare various values 
from the literature to the value derived from the dust we use in our models (see Table [J)- Our 



dust model is a combinati on of 0H5 dust ( Ossenkopf fc HenningI Il994j ) and iPoUack et al. 



(jl994 ) dust as described in Young &: Evans ( 2005 ). For 0H5 dust, we calculate the value of 
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Qa{UV) /Qe{l25fim) assuming a 10, OOOK blackbody. Qa is the stellar flux weighted average 
absorption efficiency of the dust and cTabsl-^) is the absorption cross-section for a dust grain. 
Therefore, 

QaiUV) 



Qe{l25fi'm) 



E 



X=0.15fj.m 



crabs(125/im) = 253. 



(9) 



Now we consider the various dust models in Table [T] and fix the value of jSe, the dust 



grain's emission efficiency exponent. From equation 
proffie becomes 



the form of the dust temperature 
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where f3e is fixed and the values of K{j3e) for each model are given in Table [H In Figure [1] we 
compare 0H5 dust to an analytic approximation (see equation [5]) normalized at 125/im and 
vary the value of (3e. The line with (3^ = 1-8 fits well at long wavelengths but not at shorter 
wavelengths. The opposite is true for = 1-0. Since we expect long wavelength emission 
to dominate these clouds, we adopt an analytic approximation of 0H5 dust with jSe = 1.8 
as our "Analytic Solution" in the following sections. 



2.2. Numerical Dust Temperature Calculation 



Another method of calculating the dust temperature aro und multiple sources uses the 
code Dusty which we have set up with 0H5 d ust op acities (lOssenkopf fc Henning) Il994l ) 
using the method described in Young fc Eyan2 ( 2005 ). Dusty is a one-dimensional spheri- 
cal radiative transfer code (INenkova et al.l l2000l ) . Once a dust temperature distribution is 
derived around a single source, we use it to estimate the dust temperature around multiple 
sources, using some of the assumptions in the previous method, which we explain in more 
detail below. Using Dusty to calculate the dust temperature around a young star is the 
most accurate solution; however, this program can take over 1 minute to run for low optical 
depths with yet longer run times for larger r. Therefore, we calculate the dust temperature 
profile for various combinations of luminosity, outer radius, and density profile to create a 
look-up table. (The parameters we consider are discussed in §2.2. 1[ ) We then assume that 
the dust temperature profile around a single source is of the form given in equation ([8]), i.e. 



T,(Ar„) = K, 



(Ar*i' 



l/(4+ft) 



(11) 



where K and (3 are functions of the density profile and dust properties. 
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This assumption is valid when the dust is optically thin. Although the gas and dust are 
denser closer to the central source and likely to be optically thick, we are mainly interested 
in the temperature distribution far from the central source where the physical processes we 
consider are dominant and the dust and gas are optically thin to the cloud exterior. For 
example, close to a forming star [r < 1000 AU (IMcCaughrean fc 0'delllll996l )] there may be 
a disk, which is not represented by our assumption of spherical symmetry. Therefore, it is 
not useful to model the dust temperature close to the star because it will be affected by the 
presence of a disk. 



For each set of parameters, we run Dusty and solve for a value of K and (3 (see §2.2.ip . 
Since we are only interested in the dust temperature far from the source, we fit the outer 
25% of the dust temperature profile in logT — log r space using least-squares fitting in order 
to determine the values of (3 and K in equation flTTl) . We call these values of (3 and K our 
"Fit Solution" . In order to calculate the dust temperature of a region heated by more than 
one protostar, we add up the flux at the region of interest using 



N 



47r(Ar, 



(12) 



Combining equations ( ITTi) and ( fT2i) . we derive 
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In general, we assume 
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where K and (3 are defined as the flux- weighted averages of the K and (3 values that contribute 
to the flux at point r, i.e. 



and 
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where the sums are from i = 1 to the total number of stars, A^. 

Therefore, equating (|T3l) and f|T^ and solving for Trf(r) gives. 
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Using equation (fTTl) . we obtain 



l/(4+/3) 



Ar2 ' 
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the equation we used to calculate the dust temperature in a field of sources. 

2.2.1. Parameter Space 

We assume the dust and gas are well- mixed and have the same density profile, offset 
only by the dust to gas mass ratio ijidg) described in ^ i.e. pdust = VdgPgas We assume 
Nii^/Niie = 5, which gives fj, = 2.33. The density profile of the gas is parameterized with n^, 
and a using 



We choose an inner radius of the dust distribution to be fixed at 30AU. This number 



is smaller than the average sizes of more evolved disks [i.e. iMcCaughrean fc O'delll (119961 ) 
find disk diameters of 50-1000 AU for disks in Orion], which would be considered sub-grid 
physics in a simulation of clustered star formation with moderate resolution. However, we 
find that the choice of the inner radius does not greatly affect the values of f3 and K for 
our models since we are not modeling the increase in temperature where the dust becomes 
optically thick. 

The luminosity required to reach a dust destruction temperature of 1500 K at 30AU 
using the Analytic Solution is 9.58 x lO^L©. This is outside our sample range of luminosities; 
therefore we can neglect the effects of dust destruction. We also vary the outer radius, rout, 
to obtain a range of values. We find that our values of K and /3 do not strongly depend on 
rout- 

In all cases, the stellar input spectrum is assumed to be a blackbody with temperature 
T^: = 10, OOOK. This value does not strongly infiuence the output temperature distribution 
at outer radii since the light is quickly reprocessed by the dust to longer wavelengths. 

We model the entire parameter space listed in Table [2] with two exceptions. First, we 
limit the density at the inner radius to be less than 10^°cm~^, because at higher densities 
dynamical effects may become important since the free-fall time becomes small as the density 
increases. Second, we limit the mass of the envelope to be less than ~ 1000 Mq since a larger 
envelope would likely produce a cluster of stars (assuming a star formation efficiency of 10% 
and a maximum stellar mass of 100 Mq) which would break the assumption of spherical 
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symmetry. Therefore, combinations of large ot and large Uo may not be represented in our 
models. Based on these restrictions, of the 5049 possible models in our parameter space, we 
model 3231 or 64% of them. In Figure [2] we show the luminosities and masses modeled in 
our parameter space. Figure [3] shows the relationship between the values of r, a and for 
the models in our parameter space. 

2.3. Comparison of Dust Temperature Calculation Methods 

Figures m and [5] compare the two methods of calculating the dust temperature around a 
single source. The Analytic Solution uses 0H5 dust parameters with /3e = 1.8 as described 
in §2.11 which best matches the dust properties used in the Dusty Solution in Figure [H 

Figure H] shows that the Analytic Solution captures the shape of the Dusty temperature 
profile, but at high r it overestimates the magnitude. If we only vary the luminosity (see 
Figure E]) the offset of the Analytic Solution to the Dusty Solution changes, therefore a simple 
adjustment to the Analytic Solution based on luminosity would not correct the Analytic 
Solution. The Fit Solution described in §2.21 appears to be the best fit, as shown in Figures 
Hand El 

Figure [6] shows histograms of K and /3 derived from our Fit Solution for our parameter 
space. While most models cluster around the Analytic Solution, there is a spread that is 
dependent on some of the input parameters. In Figure [7] we show how the values of K and /5 
depend on the a parameter. Models that are far from the Analytic Solution (i.e. K = 1000, 
/3 = 0) tend to have small values of a. 

In Figures [8] and [9l we show how different parameters determine the values of K and (5. 
Figure [8] shows that luminosity and the value of K are positively correlated. Although we 
have explicitly removed the luminosity dependence from equation ffTTj) . there is still some 
dependence of K on luminosity. This can be understood in terms of radiative trapping. For 
high luminosities, more photons at shorter wavelengths exist farther from the star, which 
increases the size of the region of high optical depth and therefore increases the value of 
K. Also, as a increases, K decreases. Increasing a shrinks the size of the region of high 
optical depth due to the buildup of material close to the star which absorbs the high energy 
photons. Figure M shows that as a decreases and luminosity increases, f3 moves farther 
from the Analytic Solution. The drastic drop of (3 at high luminosities can be understood by 
comparing the two models shown in Figure O In the two models, the radius at which the dust 
temperature turns up (i.e., when the material changes from optically thin to optically thick) 
moves out in radius as the luminosity increases. Since we model the value of (3 using only 
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the outer 25% of the material, we expect to be considering only the optically thin material. 
However, as the luminosity increases, the radius at which the transition from optically thick 
to optically thin material moves outward. Therefore, at higher luminosities, our calculation 
of (3 becomes influenced by the optically thick region. This is the reason (3 moves away from 
the optically thin Analytic Solution as luminosity increases. 

In order to calculate values of (3 and K not modeled in our parameter space, we inter- 
polate between known values in our look-up table . We use the met hod for interpolating in 



two or more dimensions described in Chapter 3 of iPress et al.l (119921 ). This method involves 
solving successive one-dimensional interpolations. We modified the P0LIN2 subroutine to 
interpolate in 4 dimensions. The actual method of interpolation that we used was the 
polynomial i nterp olation method over 3 known quantities from the subroutine POLINT in 



Press et al.l (Il992l ). We tested our interpolation method by calculating extra models, not 
included in our model grid, and comparing the results between the real solution and the in- 
terpolated value. The results are tabulated in Table [31 We have varied all of the parameters 
by different amounts and found that the largest percent difference from the model to the 
interpolated value is 6.0%. The largest differences occur at the extremes of the parameter- 
space at small a and large luminosity. Large errors at small values of a are likely due to the 
large scatter of values of K and at small a (see Figures [8] and [9]) . The large differences 
at large luminosities are probably due to the rapid change of the value of (3 with luminosity 
(see Figure [9]). Although K does not change as much as /?, the value of K is dependent upon 
the value of [3. 

Based on the previous discussion in this section, we use the Numerical Dust Temperature 
method described in §2.21 to calculate the dust temperature in the remainder of the paper. 
Since we are primarily interested in the dust temperature far from the luminosity source, 
we find that we can model the shape and magnitude of the dust temperature distribution 
most accurately with the Numerical Dust Temperature method. There might be some error 
in the calculation of K and (3 for models with large lu minosities or small a 's. However, 



we don't expect very small values of a to be common (iMueller et al.l |2002| ) . Also, stars 
that will eventually have high luminosities are rare and will take awhile to reach this state, 
therefore we do not expect many stars to have high luminosities during the early stages of 
star formation which we model. 



3. Gas Temperature Calculation 



After the dust temperature as a function of distance from luminosity sources is derived 
for positions near stars in a cluster using the look-up table, the gas temperature can be 
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calculated assuming gas energetics balance. We calculate the gas temperature using a gas- 
dust energetics code which includes energy tran sfer between gas and dust via collisions, 
he ating by cosm i c rays , and molecular coohng (see iDoty &: Neufeldl ( 119971 ) and the appendix 



m 



Young et al.l (120041 ) for a m ore detailed description). W e assume that the dust to gas 



mass ratio is r]dg = 4.86 x 10 ^ (IHol 
bary on is 6.09 x 10^^^ ([Young et al. 



enbach fc McKedll989l ) and the grain cross-section per 
20041 ). The cosmic ray ionization rate is 3.00 x 10~^^ 



s~^ (Ivan der Tak fc van Dishoeckl 120001) and the energy deposited per cosmic ray ionization 
is 2.00 X 10^ eV ( Goldsmith! 2001 ) in our models. We take the fraction a l abun dance of CO 
relative to H2 to be 1.0 x 10~^ from Figure 8a of iLee. Bergin. &: Evand (120041 ). We assume 
that the region we are studying is deep within a larger molecular cloud. Therefore, there is no 
interstellar radiation field impinging on the outer bounds of the cloud and the photoelectric 
effect on PAH's is not present. 

The model dependent input parameters are Tdust, local density, column density, and 
local velocity dispersion (b). Tdust is calculated with the procedure described in §2.2[ Local 
density and column density can be derived from our input density profile. The column 
density is calculated radially from the point of interest to the edge of the system. The 
edge of the system is defined as either the point at which the density is lowest or some 
fiducial value (as discussed later in §1]). The velocit y-spread para meter, b, is defined for a 
Maxwellian velocity distribution as b = (2fcT/m)^/^ (ISpitzerlll998l ) and is assumed to be 1 
km/s throughout this paper. 

Figure [TU] shows the variation of gas temperature with distance from a stellar heating 
source for two values of X(CO). Close to the source, the dust and gas temperature are coupled 
due to coUisional interactions of the dust with the gas. As the density decreases, collisions 
between the dust and gas become less frequent and the gas is able to cool via molecular 
(mainly CO) rotational transitions. Then as the density continues to drop and there is less 
CO to cool the gas, cosmic ray heating becomes the dominant heating source and the gas 
temperature increases. For a gas with less CO (X(CO)= 5 x 10~^), the cooling is not as 
efficient and the temperature is larger. As various other parameters change in our models, 
different heating and cooling terms dominate and the minimum and maximum temperatures 
vary. We discuss this in more detail in the following section (§!]). 

As seen in Figure [TOl the gas temperature falls below lOK when gas collisions are 
the dominant coolin g meth od. Our gas cooling rate calcu l ations are based on the work of 
Neufeld fc KaufmanI (119931 ) and iNeufeld. Lepp. fc Melnickl (119951 ) which only extend down 
to lOK because the H2-CO coUisional rates were not defined below lOK at the time of their 
work. We expect the cooling rate to drop drastically as the temperature approaches zero 
and we have attempted to adjust our coohng rate calculation to account for this. Our first 
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attempt to calculate the cooling rate (A) below lOK involved calculating the rate of change of 
log A with log T between 20K and lOK (cilogA/ci l ogT) and applying it to temperatures below 
lOK. This method has been used in lYoung et al.l (120041 ) successfully when dust-gas collisions, 
rather than H2-CO collisions, are the dominant mode of energy transfer at low temperatures. 
However, for our models, H2-CO collisions dominate at low temperatures and to take this 
into account we modify the method of calculating dlogA/dlogT using the Large Velocity 
Gradient model (LVG). Instead of extra polating the cooling rat e to low temperatures, we 
extrapolate the CO rate coefficients from lFlower fc Launayl (119851 ) from lOK to 5K. Then we 
calculate c/logA/rflogT between lOK and 5K using the LVG model. We apply the new values 
of (ilogA/dlogT between lOK and 5K to our gas energetics model. In Figure [TO] we compare 
the temperature derived using the old and new method of calculating the gas temperature. 
Our change increases the temperature by approximately IK at the minimum value. 



4. Gas and Dust Temperature with Multiple Sources 
4.1. Two Sources 

In order to calculate the gas temperature between two sources, we must first calculate 
the dust temperature. We do this using the formalism discussed in §2.2[ Once we have 
determined the dust temperature, we can use our energetics algorithm to calculate the gas 
temperature. Around each source we place a density profile. In order for this to be realistic, 
we choose a density, at which we have the two density profiles meet. The value of n^q sets 
the distance between the sources, i.e. smaller values of place the sources farther apart. 

Figure [11] shows the dust and gas temperature profile for increasing values of Ueq, i.e. 
smaller separations. An interesting feature of this plot is that if we only look at the region 
between the two sources, we find a large variation in gas temperature. This is due to the 
higher densities sampled as the sources move closer together. The top panel shows sources 
that are far apart and we see that the maximum gas temperature between the two sources 
is ~ 20K and the minimum gas temperature is ~ 7K. Cosmic ray heating, though relatively 
weak, can warm the material sufficiently far from luminous sources. As the sources move 
closer together, cosmic rays become less important until the temperature between the sources 
ceases to increase, whereas dust heating becomes more important and the gas is not able to 
cool efficiently and the minimum temperature between the sources rises. 
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4.2. Three Sources 

Here we calculate the gas and dust temperature distribution around three sources. The 
three sources were placed on three of the corners of a square with sides of length lOOOAU. 
The least luminous source was placed on the corner between the other two sources. The 
parameters assumed for each source are given in Table H] as well as the position of the 
sources in a 2000AU x 2000AU grid. We calculate the dust temperature using the method 
described in §2.21 and show the results in Figure [T2j 

In order to determine the gas temperature, we first calculate the density (p) and column 
density (Ncoi) at each point in the grid due to all three sources, individually. Then, at each 
point we choose the source which gives the highest value of Ncoi and use that source to 
calculate Ncoi, p, and the gas temperature. 

In order to calculate Ncoi, we integrate from the point of interest to the "edge" in the 
direction radially away from the source. We tested two methods of defining the "edge." 
Our first method, the "Length of Square Method," integrates from the point of interest to 
2000AU from the source (Figure [T3] and [T^ . The value of 2000 AU was chosen to equal the 
length of the side of the square in which we place our sources. Our second method, the "Edge 
of Square Method," integrates from the point of interest to the edge of the region studied 
(Figure (TB] and [TB]) . Using this method, the integration length depends on the direction 
of integration. In both cases, we find that the gas and dust temperature are not equal, 
unless we are in a region of high density close to a luminosity source as seen in Figures [13] 
and [151 In these high density regions, the dust and gas temperature are coupled through 
collisions. As the density decreases, the gas temperature drops rapidly, compared to the dust 
temperature, due to the ability of the gas to cool through molecular transitions. Another 
interesting feature of these plots is the detectability of Source 1, the dimmest object in the 
region, even though it is close to Source 2, the brightest source in the region. 

Two differences between the two methods of calculating N^oi are evident in Figures [T^- 
\T5[ The first is the the square shape of the contours near the edges when using the Edge of 
Square Method in Figure [T5l This is an artifact of a square region of interest. The second 
difference is the gas temperature near Source 1. It appears to be shifted in the direction of 
Source 2 for the Length of Square Method (Figure [T3l) and is shifted in the opposite direction 
for the Edge of Square Method (Figure [T^. For both methods, the dust temperature and 
gas-dust coUisional heating is the same in this region. Therefore, the molecular cooling rate 
and the column density must be different. For the Length of Square Method the calculation 
of Ncoi is dominated by Source 2 near Source 1. This is the reason the gas temperature 
contour around Source 1 is shifted toward Source 2. The value of Ncoi drops as we move 
across Source 1. This causes a shift in the gas temperature toward Source 2 even though the 
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dust temperature on either side of Source 1 is symmetric. For the Edge of Square Method, 
Source 1 is able to increase the value of N^i on the side facing away from Source 2 over the 
value calculated from Source 2. This increase in N^oi causes the gas temperature to increase 
as well. Although neither of these methods is entirely correct, we use the Length of Square 
Method (adjusted for the size of the square) in the following figures due to the square edge 
effects of the Edge of Square Method. 

In Figure [T71 we compare the gas and dust temperature far from the sources when we 
replace the three sources with a single source at the location of Source 2. We use the same 
parameters as Source 2 for our single source, however we have increased its luminosity to 
130^0. This figure illustrates the difficulty in determining the number of sources responsible 
for heating the gas and dust and highlights the need for adequate spatial resolution. 

In Figure [TSl we have zoomed out of the region of interest. Notice that the dust 
temperature steadily decreases as the distance from the central three sources increases. Yet, 
the gas temperature slowly begins to rise. This is due to the decrease of effectiveness of CO 
cooling and the increase in heating by cosmic rays as the density decreases. 

5. Conclusion 

We have presented a method for calculating the dust and gas temperature between stellar 
sources. The analytic method that we investigated for calculating the dust temperature was 
not accurate enough. Instead, our chosen method of calculating the dust temperature uses a 
simple radiative transfer code which we use to create a look-up table. Once we have derived 
the dust temperature, we are able to calculate the gas temperature by balancing various 
energy processes. We include dust-gas coUisional heating, molecular cooling, and cosmic-ray 
heating. When we have balanced the energies, we are able to derive the gas temperature. 
Other methods which set the gas temperature and dust temperature equal assume the gas 
and dust are opaque everywhere, which is not always true. In Figure dni we show the 
percentage difference between the gas and dust temperature, as well as the density, for the 
distribution of sources discussed in §4.21 These two figures show that the largest percentage 
difference between the dust and gas temperature occurs where the density is the lowest. 
Therefore, at low densities (n < lO^cm^^), assuming equal dust and gas temperatures is not 
appropriate. 



We plan to use the method discussed in this paper to model a re gion of clustered sta r 
formation with the three-dimensional hydrodynamics code discussed in lMartel et al.l (120061 ). 
Our method of calculating the gas and dust temperature distribution in a field of young stars 
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will enable us and others to more accurately model clustered star formation observationally 
and in future simulations. 
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Table 1. Dust Parameters 



Dust Type 


Qa.{UV) 

Qe(125/jm) 


K{1) 


K(1.8) 


K{2) 


OH5 (our dust model) 


253 


539.3 


351.1 


319.8 


Hildebrand (1983) 


4000 


936.6 


565.1 


506.6 


Makinen et al. (1985) 


790 


677.1 


427.2 


386.7 



Table 2. Dust Model Parameters 



Parameter 


Lower Limit 


Upper Limit 






log (L/Lq) 


-2 


6 


0.5 


17 


log (no/lcm~^) 


2.5 


7.5 


0.5 


11 


a 





4 


0.5 


9 


log (rout/1 pc) 


-1 





0.5 


3 



^Spacing of parameters 
'^Number of parameters 



Table 3. Dust Model Interpolation Results 



L(Lq) 


rio (cm ^) 


a 


Tout (pc) 




~^poly — inter p 


]^%-dtff 




l^poly — inter p 




1.8 X 10-2 


5.62 X 10^ 


1.75 


0.177 


302.1 


300.8 


0.4 


1.82 


1.82 


0.0 


1.0 X 10° 


5.00 X 10^ 


2.75 


0.194 


214.8 


213.9 


0.4 


1.81 


1.81 


0.0 


1.0 X 10° 


5.00 X 10^ 


1.75 


0.194 


295.6 


295.1 


0.2 


1.78 


1.78 


0.0 


1.0 X 10° 


5.00 X 102 


0.80 


0.194 


322.4 


310.0 


3.8 


1.68 


1.78 


6.0 


1.0 X 10" 


5.50 X 10^ 


2.00 


0.194 


161.5 


157.6 


2.4 


1.82 


1.81 


0.5 


1.5 X 10^ 


5.00 X 10^ 


2.00 


0.194 


355.8 


356.5 


0.2 


1.37 


1.36 


0.7 



Table 4. Source Parameters 



Source 


Luminosity [Lq) 


no(cm ) 


a 


rout (pc) 


K 


/3 


X (AU) 


y (AU) 


1 


1 


10^ 


2 


0.1 


302.60 


1.7875 


500 


1500 


2 


100 


10^ 


2 


0.1 


215.89 


1.7623 


1500 


1500 


3 


10 


10* 


2 


0.1 


263.22 


1.7817 


500 


500 
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Fig. 1. — 0H5 dust properties. Dashed line shows the variation of cross-section with wave- 
length for 0H5 dust. The solid lines show the different values of jSe normalized at 125/im 
that we consider in Table [H 
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Fig. 2. — Masses and luminosities of models in our parameter space. Figure on the left 
shows a histogram of the masses of the clouds surrounding a source. On the right, the 
relationship between the mass and luminosity of all our models is plotted. For every density 
distribution in our sample, we have models which correspond to all of the luminosity values 
in our parameter space. 
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Fig. 3. — Relationship between a, t, and Hq. The figure on the left shows the relationship 
between a and Ug for all our models. Some models at low a with high values of Ug are not 
in our sample due to the maximum mass criterion. Models missing in the top right corner 
of the figure at high Up and high a are missing due to the criterion which sets the maximum 
density at the inner radius. On the right, we show how r varies with a, Uo, and Tout for our 
models. At a fixed value of Uo, as a increases (and a > 1), t increases as well due to the 
sharp density increase in the profile. For lower values of a (a < 1), r begins to increase 
again, but this increase depends on amount of material included in the profile at the edge, 
i.e. the value of rout- 
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Fig. 4. — Comparison between temperature distributions for two different models with a low 
and high fiducial density. For both models, L = ILq, a = 2, and r^^ut = 1 pc- The figure 
on the left has log Uo = 2.5 and r = 2.894 x 10^"^. The fit parameters are K = 300.79 and 
(3 = 1.822. The figure on the right has log no = 5.5 and r = 2.894 x 10"^ with K = 169.5 
and P — 1.799 . Although both of these models are optically thin, as assumed in the Analytic 
Solution, it is clear that it is not a good description of the dust temperature. If we were to 
change the value of K assumed in the Analytic Solution, we might be able to fit the exact 
solution provided by Dusty for one model, but it would then not fit for another model. 



-22- 




io> 



10* 



I 

log r (AU) 



10* 



10> 



■ ■■■I 

10« 



I 1_1_L. 

log r (AU) 



10* 



Fig. 5. — Comparison between temperature distributions for two different models with a low 
and high luminosity. For both models, a = 3, and Tout = O.f pc, log Uo = 2.5. The figure 
on the left has log L/Lq = —1. The fit parameters are K = 276.23 and (3 = f.802 The 
figure on the right has log L/Lq = 3 with K — 297.5 and /3 = 1.658. These models show 
how increasing the luminosity increases the overall dust temperature as well as making the 
optically thick region near the center extend farther out. 




Fig. 6. — Figures shows the range of K and j3 for the chosen parameters in Table O The 
vertical line marks the values of K and (3^ in the Analytic Solution. 




Fig. 7. — The range of K and (3 for the chosen parameters in Table [2] as a function of a. The 
vertical line marks the values of K and /5e in the Analytic Solution. Histograms for a = 
have been multiplied by 2 for clarity. 
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Fig. 8. — The horizontal hne marks the value of K in the Analytic Solution. Individual 
models are plotted as circles of various sizes. The size of the circle indicates the value of 
the a parameter as noted in the top of the plot. The nine separate plots each show K a.s a. 
function of Uo for nine different luminosity regimes. In the top-left box, log L/Lq is —2 or 
— 1.5, as indicated at the top- right corner in the box. The bottom- right box shows models 
with the highest luminosities. 
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Fig. 9. — The horizontal hne marks the value of (3^ in the Analytic Solution. Same plot 
details as Figure [8] except /5 is plotted rather than K. 
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Fig. 10. — Gas temperature distribution and comparison of old and new cooling rates. Solid 
line shows Fit Solution to dust temperature. The model parameters are log {L/Lq — 1), log 
(^^out/lpc) = —0.5, log n„ = 4.5, and a = 2.5. Dashed and dotted line show gas temperature 
for X(C0) = 1 X 10"^ and X(C0)=5 x 10~^, respectively. Old and new cooling rates are 
shown in the thin and thick lines, respectively. 
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Fig. 11. — Gas temperature distribution between two sources. Plots show dust (thick hue) 
and gas (thin hue) temperature as a function of distance between two stellar sources. From 
top to bottom, distance between sources is decreasing, such that the gas density surrounding 
the two sources agrees with the value, Ueq quoted on the left. The source on the left (1) 
has the parameters L — IL©, Ug = 10^, and a — 2. The right source (2) has L — lO^Lo, 
Ug — 10^, and a — 2. The three horizontal plots for the different values of Ueg show the same 
data from different perspectives. The large plot on the left is plotted on a linear scale from 
Source 1. The two smaller plots on the right are plotted on a logarithmic scale, from the 
perspective of Source 1 (middle) and Source 2 (right). This is done to show structure close 
to the individual sources. 




Fig. 12. — Surface plots of dust temperature. Sources are labeled according to the parameters 
listed in Table HI The first figure shows the full range of dust temperatures and the second 
figure shows a smaller range of dust temperatures. 
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Fig. 13. — Contour plots of dust and gas temperature using the Length of Square Method. 
Sources are labeled according to the parameters listed in Table ID 
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Fig. 14. — Surface plots of gas temperature using the Length of Square Method. The same 
data is shown in both plots with a smaller range of gas temperature shown in the second 
plot. 
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Fig. 15. — Contour plots of dust and gas temperature using the Edge of Square Method. 
Sources are labeled according to the parameters listed in Table HI 
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Fig. 16. — Surface plots of gas temperature using the Edge of Square Method. The same 
data is shown in both plots with a smaller range of gas temperature shown in the second 
plot. 




Fig. 17. — Contour plots of dust and gas temperature with three sources and one source. 
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Dust Temperature 




Fig. 18. — Surface plots of dust and gas temperature. Top plots show dust temperature. 
Bottom plots show gas temperature. Plot on top left show complete sample of dust tem- 
peratures. Top right plot shows a smaller range in dust temperatures. Bottom right and 
bottom left show same data from two different orientations. Area sampled is 50000 x 50000 
AU. The Length of Square Method is used here. 
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Fig. 19. — Percentage difference between gas and dust temperature and density. The same 
data and method as seen in Figures [T^ - [TH First figure shows percentage difference of dust 
and gas. Second figure shows density. 



